/************************************************************************
*************************************************************************
*************************************************************************

gendatsyn.do

Code to generate a hypothetical dataset, with financial constraints and/or
variation in input prices. This version uses syntax rather than an ordered
list of arguments

2 March 2016: 		Created




The production function is
X(K,L)*M^pm

OR

X(K,L)*(pm1*M - M^2)

Arguments:
klfunc mfunc constraintmult rho pk pl pm N T sigeps sigom pricemsd

klfunc			-	a command that uses K and L to generate X(K,L) and stores it in a variable X_kl
mfunc 			- 	if equal to "quad", use the quadratic form; if equal to "exp", use the exponential form; otherwise use isoelastic
constraintmult	-	the multiplier on "wealth" in the liquidity constraint
rho				-	persistence of productivity
N				-	number of firms
T				-	number of years in the dataset
pm pm2 			- 	The cobb douglas coefficient, or if mfunc=quad the coefficient on the squared term
sigeps			-	standard deviation of the unanticipated shock
sigeta			-	standard deviation of the innovation to productivity
pricemsd		-	standard deviation of the price of flexible inputs
meanprice		-	average price (Default: 1)
expshift		-	Shifter in choice of intermediates: Y = exp(om+skh) X_kl (1 - exp( -theta*(M-expshift) ) )
					[default: value chosen to ensure no corner solutions]





There are also the parameters of the stochastic processes for capital, labor,
and wealth
	
	
	/*
	These constants are calibrated from Industry 311 in the Chilean ENIA
	*/
	
	Capital:
	k 			= `kc' + `kk1'*l.k + `ke1'*l.eta  + `ke2'*l2.eta 	if t > 2
	local kk1 = 1.003592
	local ke1 = .0359078
	local ke2 = .1957049
	local kc  = -.012955
	
	Labor:
	l 			= `lc' + `ll1'*l.l + `le1'*l.eta  + `le2'*l2.eta	if t > 2
	local ll1 = .9519013
	local le1 = .4676349
	local le2 = .1845552
	local lc  = .188374
	
	
	Mean and Std. Dev:
	gen k 			= 2.115286*rnormal() 	+ 6.946116 	
	gen l 			= 0.9518106*rnormal() 	+ 3.659235	

	
	
		/*
	These constants are calibrated from Industry 311 in the Chilean ENIA
	(I treat log[stassets] as S )
	*/

	Savings:
	gen savings 		= `savmean' + rnormal()* `savsd'
	replace savings 	= `savc' + `savr1'*l.savings + `savr2'*l2.savings + rnormal()* `savres' if t > 2
	
	
	local savr1			= .4920887
	local savr2 		= .3205638
	local savc 			= 1.940491
	
	local savmean		= 7.417812			//Overall mean
	local savsd 		= 2.438074 			//Total standard deviation
	
	local savres 		= 1.354485 			//Standard deviation of the residual
	
	
	Unobservable:
	(If any parameter is not specified, we assume it is identical to the value
	for "savings")
	
	foreach piece in r1 r2 c mean sd res {
		if mi("`unobs`piece''") {
			local unobs`piece' = sav`piece'
		}
	}
	
	mmerr: Fraction of standard deviation that consists of iid measurement error;
			This is added only to k,l,m (as y already has epsilon)
			
	constraint: The function used to determine the constraint, which will
					be set equal to 
					`constraintmult' * collateral^(1-`constraintmix') * `constraintvar'^`constraintmix'

*************************************************************************
*************************************************************************
************************************************************************/


cap program drop gendatsyn
program gendatsyn
	noisily: di "in prog"
	#delimit;
	syntax  =/exp , klfunc(string) rho(real) N(integer) T(integer) pm(real) sigeps(real) sigeta(real)
				[	pricemsd(real 0) nlprod(string) fracobs(real .5)
					kk1(real 1.003592) kk2(real 0) ke1(real .0359078) ke2(real .1957049) kc(real -.012955) 
					ll1(real .9519013) ll2(real 0) le1(real .4676349) le2(real .1845552) lc(real .188374)
					ksd(real 2.115286) kmean(real 6.946116) lsd(real 0.9518106) lmean(real 3.659235	)
					savmean(real 7.417812) savr1(real .4920887) savr2(real .3205638) 
					savc(real 1.940491) savsd(real 2.438074) savres(real 1.354485)
					unobsmean(real -1) unobsr1(real -1) unobsr2(real -1) 
					unobsc(real -1) unobssd(real -1) unobsres(real -1) 
					mmerr(real 0) mmerrcor(real 0) constraintvar(string) constraintmix(real 0) 
					mfunc(string) meanprice(real 1) expshift(real -1)
				];
	#delimit cr

	noisily: di "post syntax"
	//The initial expression is the multiplier on the constraint
	local constraintmult = `exp'
	
	//Stata won't let us have capital letters in macros passed; let's set them now
	local N = `n'
	local T = `t'
	
	

	//////////////////////
	//////Set up firms and productivity
	
	//Variance of productivity
	//local sigom = sqrt( (1-`rho')^2 / (1-`rho'^2) ) * `sigeta'
	local sigom = sqrt( 1 / (1-`rho'^2) ) * `sigeta'
	
	//Initial state
	clear
	set obs `N'
	gen id 			= _n
	gen omega 		= rnormal()*`sigom'
	gen k 			= `ksd'*rnormal() 	+ `kmean' 	
	gen l 			= `lsd'*rnormal() 	+ `lmean'	
	gen logpricem 	= log( `meanprice' )
	
	//Hypothetical fixed constraint
	gen confe		= exp( rnormal() )
	
	//Add time dimension
	expand `=`T'+2'
	bysort id: gen t = _n
	
	xtset id t
	
	//////////////////////
	//////Create the flow of productivity and factors
	
	//Random shock and its mean
	if !mi("`sigeps'") {
		gen eps 		= rnormal()*`sigeps'
	}
	else {
		gen eps 		= rnormal()
	}
	gen tmp 		= exp(eps)
	quietly: sum tmp
	gen meaneps 	= r(mean)
	drop tmp
	
	di "`pricemsd'"
	
	//Generate input price shocks
	if !mi("`pricemsd'") {
		gen priceeps 		= rnormal()*`pricemsd'
	}
	else {
		gen priceeps 		= rnormal()*0
	}
	
	//Generate future capital, labor, and productivity
	//gen rainshock 		= rnormal()*`sigeta'											if t > 1
	//gen eta 			= .5*rainshock + rnormal()*`sigeta'
	
	gen eta 				= rnormal()*`sigeta'
	
	if mi("`nlprod'") {
		//replace omega 		= `rho'*l.omega + (1-`rho')*eta 							if t > 1
		replace omega 		= `rho'*l.omega + eta 										if t > 1
	}
	else {
		replace omega 		= `nlprod' 	 + eta 											if t > 1
	}
	
	
	
	
	
	
	
	replace k 			= `kc' + `kk1'*l.k + `kk2'*l2.k + `ke1'*l.eta  + `ke2'*l2.eta 	if t > 2
	replace l 			= `lc' + `ll1'*l.l + `ll2'*l2.l + `le1'*l.eta  + `le2'*l2.eta	if t > 2
	
	
	//Persistence of prices (calibrated from data)
	local persistp .9274344
	
	local pricescale = sqrt(1-`persistp'^2)											//Rescale error to keep the variance stationary
	replace logpricem	= `persistp'*l.logpricem + `pricescale'*priceeps	if t > 1
	
	gen K 				= exp(k)
	gen L 				= exp(l)
	gen A				= exp(omega)
	
	

	gen savings 		= `savmean' + rnormal()* `savsd'
	replace savings 	= `savc' + `savr1'*l.savings + `savr2'*l2.savings + rnormal()* `savres' if t > 2
	
	/*
	(If any parameter of the stochastic process for the "unobservable" is not 
	specified, we assume it is identical to the value for "savings")
	*/
	foreach piece in r1 r2 c mean sd res {
		if `unobs`piece'' < 0 {
			local unobs`piece' = `sav`piece''
		}
	}
	
	gen unobs 			= `unobsmean' + rnormal()* `unobssd'
	replace unobs 		= `unobsc' + `unobsr1'*l.unobs + `unobsr2'*l2.unobs + rnormal()* `unobsres' if t > 2
	
	
	gen wealth 			= exp( sqrt(`fracobs')*savings + sqrt((1-`fracobs'))*unobs )
	
	//Rescale by means
	quietly: sum K
	local scK = r(mean)
	
	quietly: sum wealth
	local scW = r(mean)
	
	
	gen collateral		= (K/`scK')^.5 *  (wealth/`scW')^.5
	
	//If they haven't specified a constraint, just assume it's the collateral
	if mi("`constraintvar'") {
		gen con				= `constraintmult' * collateral
	}
	else {
		
		//If there's an operator, generate a new variable to go with it
		if regexm("`constraintvar'","\.") {
			tempvar tmp
			gen `tmp' = `constraintvar'
			local constraintvar `tmp'
		}
		
		quietly: sum `constraintvar'
		local scC = r(mean)
		gen con				= `constraintmult' * collateral^(1-`constraintmix') * (`constraintvar'/`scC')^`constraintmix'
	}
	
	//Price of intermediates
	gen pricem			= exp(logpricem)
	
	//Contribution of capital and labor
	`klfunc'
	
	//////////////////////
	//////Choose flexible input and compute output
	
	
	if "`mfunc'" == "quad" {
		gen Mopt		= max( 1/2 * ( `pm' - pricem / (A * meaneps * X_kl)  ) , 0 )
		gen M 			= min(  Mopt  , con / pricem )   ///
								
		gen Y 			= exp(eps + omega) * X_kl *(`pm'*M - M^2)					
	}
	else if "`mfunc'" == "exp" {
		gen Mopt		= 1/`pm' * ( omega + log(meaneps) + log(X_kl) + log(`pm') - log(pricem) )
		
		//If expshift is unspecified, we find a value that ensures strictly positive Mopt
		if `expshift' < 0 {
			sum Mopt
			local expshift = - min( r(min), 0 ) * 1.1
		}
		
		//Shift, then apply zero lower bound
		replace Mopt 	= max( Mopt + `expshift' , `expshift' )
			
		gen M 			= min(  Mopt  , con / pricem )   ///
								
		gen Y 			= exp(eps + omega) * X_kl * (1 - exp( -`pm'*( M - `expshift' ) ) )
	}
	else {
		gen Mopt 		= ( `pm' * A * meaneps * X_kl / pricem ) ^ ( 1 / (1-`pm') )
		gen M 			= min( Mopt , con / pricem )
		gen Y 			= exp(eps + omega) * X_kl * M^`pm'
	}
	
	foreach var of varlist M Y Mopt {
		local tmp = lower("`var'")
		gen `tmp' = log(`var')
	}
	

	
	//Is the firm constrained?
	gen constrained = abs(con / pricem - M) < .0001
	
	
	if `mmerr' > 0 {
		foreach var of varlist k l m {
			gen `var'_true 	= `var'
			
			tempvar err
			
			sum `var'_true
			gen 		`err' = `r(sd)'*`mmerr'*rnormal() 
			replace		`err' = `mmerrcor'*l.`err' + (1-`mmerrcor') * `err' if t > 1
			
			
			replace `var' 	= `var' + `err'
			
			local lev = upper("`var'")
			replace `lev' = exp(`var')
		}
	}
	
	
	order y k l m, last
	
	//gen trueshare = pricem*M/Y
	
	//Throw away the initial periods
	drop if t <= 2

end



